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We theoretically study the stability conditions of the ferroelectric ice of the Cmc2i structure, 
which has been considered, for decades, one of the most promising candidates of the low temperature 
proton-ordered phase of pure ice Ih. It turned out that the Cmc2i structure is stable only with a 
certain amount of dopant and the true proton-ordered phase of pure ice Ih remains to be found at 
lower temperature. Implication for spin ice is mentioned. 
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Water is a common molecule in the universe, found 
on the earth and other solar/extrasolar planets[l, 2]. 
The solid form of water, ice, is known to have an ex- 
tremely rich phase diagram despite its simple molecular 
structure [3]. The complexity comes from configurations 
of its hydrogen bond network. Among its many phases 
ice Ih is the most abundant on earth. It is character- 
ized by hexagonal symmetry and disordered tetrahedral 
hydrogen bonds which satisfy the ice rules[4]. The resid- 
ual entropy 5*0 of 3.5 (J/mol K) due to the disorder is 
observed when ice Ih is cooled toward the absolute zero. 
At a glance this may seem to contradict the third law of 
thermodynamics, which states that entropy will approach 
zero as temperature approaches absolute zero, which is 
why the low temperature proton-ordered phase of ice Ih 
is long sought [5] . The existence of proton-ordered ice in 
space [6, 7] and its role in planet formation [8, 9] have 
also been discussed by several authors. It has been exper- 
imentally found that, when doped with salt such as KOH 
(or placed in an electric field), ice Ih transforms to "ice 
XI" below 72 K [6-14], which is a proton-ordered, ferro- 
electric crystal with space group Cmc2i. The traditional 
view is that the pure Cmc2i structure is thermodynam- 
ically stable at low temperature and the role of a dopant 
is that of a catalyst In this Letter, we examine the sta- 
bility conditions of the Cmc2i structure, which has been 
considered for decades one of the most promising candi- 
dates for a low temperature proton-ordered phase of pure 
ice Ih[6-24]. A simple model based on first principles cal- 
culation suggests that the Cmc2i structure is stable only 
with a certain amount of dopant, and the true proton- 
ordered phase of pure ice Ih remains to be found at lower 
temperature. 

The proton-ordered phase of ice Ih has also been the- 
oretically studied and a number of works published [15- 
24]. Among them Kuo and Singer [16] and Hirsch and 
Ojamae[17] identified crystallographically inequivalent 16 
proton-ordered structures for an orthorhombic unit cell 
with eight water molecules (Fig. 1). This classification 
enabled systematic comparison of hydrogen bond config- 
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FIG. 1: The figure shows the Cmc2i hydrogen bond config- 
uration, one of the 16 proton-ordered structures with eight 
water molecules in orthorhombic unit cell[17]. 



urations. Following this classification we have calculated 
the Kohn-Sham total energy E^s^ permanent polariza- 
tion P07 and dielectric constant e of the 16 structures (Ta- 
ble 1) by using the first principles electronic structure cal- 
culation (see below for details). The phonon contribution 
to the energy is neglected. The calculated energies agree 
with Hirsch's results. Structure 1 (space group Cmc2i) 
is ferroelectric and the most stable, which agrees with 
experimental observations. Structure 2 (space group 
Pna2i) is an antiferroelectric crystal, which Davidson 
and Morokuma[18] suggested early on as the most sta- 
ble structure. However, these results were calculated for 
an infinite crystal under periodic boundary conditions, 
thus effects of the macroscopic electric field £ due to the 
surface charge were neglected [25] while the local interac- 
tions between water molecules were correctly taken into 
account. 

In order to evaluate the effects of the macroscopic elec- 
tric field in polar crystallites, let us introduce a model 
with a cubic crystallite of size L ^ Q}^^ sandwiched by 
two hypothetical 'electrodes' with charge q = q^oi + Qimp 
and —q^ respectively (Fig. 2) where Vt is the volume of the 
unit cell, qimp is the doped charge, and q^oi is the surface 
polarization charge due to the polarization P = Pq + X^ 
where the first term is due to the permanent dipole mo- 
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TABLE I: Kohn-Sham total energy Eks, polarization Pq and 
the macroscopic electrostatic energy {27tPq /e)Q (in Hartree 
atomic units) of the crystallographically inequivalent 16 
proton-ordered structures of ice Ih for an orthorhombic unit 
cell with eight water molecules are listed where e = 1.8 is the 
dielectric constant. 



+ + + + + 



Up 



Ice 



FIG. 2: Model of ice crystallite: A cubic crystallite of size 
L is sandwiched by two hypothetical electrodes^ that repre- 
sent the surface charges due to the polarization P and due to 
concentration of dopant. 



ment of the oriented water molecules and the second term 
is due to the induced dipole moment. We use Gaussian 
units in this Letter. In the case of pure ice crystallite, 
where we have an open circuit boundary condition with 
D = E-\-4:7rP = or qimp = 0, the macroscopic field called 
depolarization fields £ = — (47r/e)Po5 is caused by the sur- 
face polarization charge qpoi = PoL'^/e^ where e = l+47rx 
is the dielectric constant. It was evaluated by the den- 
sity functional linear response theory that e = 1.8 for our 
systems. 

If a doped charge exists in the crystallite, it feels the 
force due to the macroscopic electric field £ and accumu- 
lates at the oppositely charged surface. Then the surface 



polarization charge is screened by the doped charge, 

Q = Qpol + Qimp = PL — pimpL . (1) 

where pimp is the average number density of the doped 
charge in the crystallite. The minimum density pmin 
is defined as the number density of the doped charge 
that screens the surface polarization charge perfectly, i.e., 
^=0 or (7 = 0, 



Po/L. 



(2) 



Dopant exceeding pmin has no effect on stability in this 
model. Let us define the dimensionless electric field x by 
S = -(47r/e)Pox or x = 1 + Qimp/qo where qo = PqP^. 
The system with x = 1 represents the pure system 
and that with x = represents the fully screened sys- 
tem. Since f = for the fully screened system it also 
corresponds to the system calculated with the conven- 
tional density functional calculation method with peri- 
odic boundary conditions, whose total energy Etot{x = 0) 
is given by the Kohn-Sham total energy Eks- In order 
to evaluate the total energy of pure system Etot{x = 1) 
we calculate the work W necessary for moving the doped 
charge q^mp across the crystallite against the macroscopic 
electric field £. Let us define infinitesimal work dW nec- 
essary for moving infinitesimal charge dqimp = qo ■ dx 
across the crystallite against the macroscopic electric 
field as 



47r 
dW = {dqimp)£L = — P^L^xdx 



(3) 



Then the total energy per unit cell at dimensionless sur- 
face charge x becomes 
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Here we find the main result of this paper that the Cmc2i 
structure of pure ice is unstable due to the electrostatic 
energy represented by the second term of eq.(4) and 
the ground state of pure ice should be non-polar, while 
the Cmc2i structure of doped ice is stable because the 
dopant acts as a stabilizer that eliminates the electro- 
static energy. All previous DFT calculations provided the 
total energy of doped ice Etotix = 0) = EkSi which is 
lower than the total energy of the pure ice Etot{x = 1) by 
the electrostatic energy. The electrostatic energy for the 
crystallite (27r/e)Po L^ is equal to the electrostatic energy 
{l/2)qQ/C of a parallel plate capacitor with capacitance 
C = (e/47r)I/. Note that the electrostatic energy is about 
hundred times larger than the variation in Eks of var- 
ious hydrogen bond configurations (see Table 1). More 
precise and sophisticated treatment of finite electric field 
calculation in the context of the self-consistent density 
functional theory is found in the recent article [25]. 



This electrostatic energy makes pure monodomain 
crystalline of typical ferroelectric materials such as 
BaTiOa unstable. According to the standard explanation 
of domain formation, the electrostatic energy of these 
materials is lowered by forming domains of polarization 
where the polarization in the half of the domains is re- 
versed to reduce the effective polarization Pg// at the 
cost of domain wall formation energy [26]. Domain size is 
determined by the balance of the gain in the electrostatic 
energy and the loss in domain formation energy. As the 
result the stable state of the pure ferroelectric crystalline 
breaks into small domains with alternating polarization 
instead of being one monodomain of crystal structure. 

In our case of pure ice crystallite, the electrostatic en- 
ergy is completely eliminated by changing the hydrogen 
bond configuration from polar to non-polar while leaving 
the oxygen lattice intact. As the result the most sta- 
ble state will be one non-polar monodomain structure. 
This hydrogen bond reconfiguration, however, can be re- 
garded as an extreme case of domain formation. For 
example, ant i- ferroelectric structure Pna2i can be con- 
sidered as ferroelectric domains with alternating direc- 
tions in molecular scale. The domain size can become so 
small because the cost of domain formation, or the en- 
ergy difference among hydrogen bond configurations, is 
extremely smaller than the gain in electrostatic energy 
(See Table 1). 
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FIG. 3: Kohn-Sham total energy of the n-th Structure 
compared to the energy of the Cmc2i structure at x = 0, 
Etot{n, x) — £^tot(l, X = 0) (in Hartree atomic units) as a func- 
tion of the dimensionless electric field x = 8/{4:7TPo/e) where 
e = 1.8 is the dielectric constant. 

Fig. 3 shows the total energy Etot{x) as a function of 
dimensionless surface charge x. The Cmc2i structure re- 
mains the most stable only up to x = 0.08. The total 
energies of the four non-polar structures (No. 2, No. 6, 
No. 7, No. 14) do not change with doping within this 
model. Therefore the most promising candidate of the 
proton-ordered form of pure ice Ih is the Structure 14 
(space group P21) which has the lowest total energy 
among the four non-polar structures. Non-polar config- 
urations in the orthorhombic unit cell containing up to 



64 molecules were explored by randomly generating hy- 
drogen bonds satisfying ice rules and then optimizing the 
geometry with density functional calculation [27]. Pre- 
liminary results do not show any non-polar configuration 
which has lower energy than the Structure 14. Random 
networks of hydrogen bond satisfying the ice rules appear 
about 200K above the energy of the Cmc2i structure. 
The phase transition temperature of doped ice, 72K, is 
much lower than 200K probably because of huge number 
of random hydrogen bond networks. Since the total en- 
ergy of the Structure 14 is lOOK higher than the Cmc2i 
structure the transition temperature of pure ice is esti- 
mated to be around 36K. 

In the case of fully doped ice, the electrostatic en- 
ergy in eq.(4) is completely eliminated by the dopant 
and the stable state of the system becomes the Cmc2i 
structure. The specimen may consist of randomly ori- 
ented crystallites (or domains) with the size constrained 
by eq.( 2). The size of crystallite is not yet experimen- 
tally well determined but it is estimated from eq.(2) to 
be as large as L = Ijim for a typical dopant density of 
p = 0.001(7710///). Even in such a case the above argu- 
ment for an isolated crystallite is considered to be valid. 
The effects of the electric fields originating from other 
domains will be negligible because the surface charge is 
screened and the crystallites orient randomly. 

Even the purest water does not consist only of H2O 
molecules but also contains hydronium ions (HsO^) and 
hydroxide ions {0H~) due to autoionization of the wa- 
ter molecules. In order to see if these ions taken into the 
ice can stabilize the Cm2i structure, let us calculate the 
minimum size of domain L^i^ = Po/^P to be stabilized 
by autoionization. Assuming ions of density 10~^ (mol/Q 
in water at standard condition are incorporated into the 
crystallite, the minimum domain size Pmin is estimated 
to be as large as 1 (cm) . Nucleation of domains of such 
large size at once would be difficult. Further, the dissoci- 
ation constant in ice is orders of magnitude smaller than 
that in water [28]. 

The details of the numerical calculations are as follow. 
The density functional electronic structure calculations 
were performed with ABINIT codes [29] based on the 
plane wave basis set, norm conserving pseudopotential, 
and GGA density functional according to Perdew, Burke 
and Ernzerhof (FEE) [30] . The Brillouin zones were sam- 
pled with the Monkhorst-Pack k-points [31] with 6x3x3 
mesh. The cut-off energy of plane wave basis was set to 
be 50 (Hartree). The positions of atoms were fully op- 
timized so as to minimize the total energy. The initial 
atomic configurations were adopted from Hirsch's table 
[17]. The unit cell was not optimized because the change 
in stress tensor due to the hydrogen bond configurations 
turned out to be very small. The dielectric constant was 
calculated by using the density functional linear response 
theory [32]. The permanent polarization was calculated 
following the Berry phase theory [33]. 

In summary we theoretically studied the stability con- 
ditions for the ferroelectric ice of the Cmc2i structure. It 



turned out that the Cmc2i structure is stable only with 
a certain amount of dopant. The true proton-ordered 
phase of pure ice Ih should be non-polar and remains to 
be found at lower temperature. We proposed the Struc- 
ture 14 (space group P21) as a promising candidate. In 
the formation of ferroelectric ice (Pq 7^ 0), dopant acts 
not only as a catalyst but also as a stabilizer eliminat- 
ing the second term of eq.(4) by screening the surface 
polarization charge. Contrarily, a dopant acts only as a 
catalyst in the formation of non-polar proton-ordered ices 
(Po = 0) [34] because the macroscopic electrostatic en- 
ergy is zero even without dopant. The recent experimen- 
tal discovery of antiferroelectric ice XV [35] in spite of the 
theoretical prediction of ferroelectric structure [36, 37] 
might be relevant to our model. Since the surface po- 
larization charge is screened by doped charges or the hy- 
drogen bond network is reconstructed to non-polar struc- 
tures, the astronomical implications of the strong electric 



field produced by ferroelectric ice [8, 9] sound unlikely. 

It has long since been clear that there is a similarity 
between water ice and spin ice [38]. In spin ice, magnetic 
monopoles[39, 40] interacting with an external magnetic 
field [41, 42] have been observed in neutron scattering ex- 
periments. In close analogy to the ferroelectric water ice, 
the effect of surface magnetic charge may also be impor- 
tant for crystallites of ferromagnetic spin ice. 
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Crystal Structure Prediction for Hydrogen Bond Materials 

We searched the hydrogen bond configurations in the 1x1x1, 2x1x1, 2x2x1 and 2x2x2 supercells of the 8 molecule 
orthorhombic unit cell [1] by using crystal structure prediction program for hydrogen bond materials. We adopted 
multi-energy scale strategy by using empirical potential to select low energy configurations and first principles calcu- 
lation to distinguish small energy differences among them. The initial configuration was generated by choosing the 
orientation of each water molecule randomly from possible six orientations, then the total energy was minimized in 
terms of water orientation by using Metropolis method for simulated annealing and PSC/E model potential [2] for 
water-water interaction. The resulting configuration was accepted if the polarization is zero and the total energy 
is less than cut-off energy to select hydrogen bond configurations that satisfy the ice-rules. We collected 100 such 
low energy configurations for each supercell. The crystallographic equivalence of the accepted configurations was not 
checked. Then the atomic positions of the accepted configurations were further optimized by using density functional 
total energy calculation program VASP[3] and the distribution of the total energy was derived. (Fig. SI). The en- 
ergy distributions were not sampled perfectly uniform but we believe they provide useful information on the energy 
distribution of hydrogen bond configurations. 
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FIG. 1: Histogram of non-polar hydrogen bond configurations with respect to the total energy of the n-th Structure compared 
to the energy of the Cmc2i structure at x = 0, Etot{ri, x = 0) — Etot{^, x = 0) (in Hartree atomic units). 
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structure No. 14 

The Structure No. 14 has the cell dimensions of a = 7.7808, b = 7.33581, c = 4.49225 in the unit of Angstrom and 
belongs to crystallographic space group P21. The fractional coordinates are listed in Table S I. 



atom X y 



HI 0.91530 0.19583 0.24980 

H2 0.79656 0.02001 0.25135 

H3 0.52649 0.51741 0.06879 

H4 0.70263 0.51818 0.24625 

H5 0.02160 0.98234 0.92673 

H6 0.02213 0.98236 0.57414 

H7 0.41772 0.30443 0.74775 

H8 0.47740 0.47955 0.56842 

09 0.91800 0.06071 0.25039 

010 0.58290 0.56414 0.25143 
Oil 0.08451 0.93553 0.75071 
012 0.41707 0.43952 0.75151 

S I: Fractional coordinates of Structure No. 14. 
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